function plotFactors(setup,output,AVarTheta1)

regime_1  = real(setup.econAct)>0;
regime_2  = real(setup.econAct)<0;

figure('Name','Annualized Factor plot','NumberTitle','off');
if setup.nn+setup.nm <= 2
    numRow = 2;
    numCol = 2;
else
    numRow = 3;
    numCol = 2;
end
dates = datenum(setup.timeIndex);
idx = 0;
for i=1:setup.nm
    idx = idx + 1;
    subplot(numRow,numCol,i)    
    hold on
    timeIndex = str2num(datestr(dates,10))+str2num(datestr(dates,5))/12;
    seriePlot = setup.macros(:,i);
    minY      = min(seriePlot);
    maxY      = max(seriePlot)*1.2;
    X  = [timeIndex(1,1):1/12:timeIndex(end,1),fliplr(timeIndex(1,1):1/12:timeIndex(end,1))];
    % For shading above zero
    Y1 = [maxY*regime_2' 0*regime_2'];
    fill(X,Y1,[0.8 0.8 0.8]);
    % For shading below zero
    Y2 = [min(0,minY*regime_2') 0*regime_2'];
    fill(X,Y2,[0.8 0.8 0.8]);
    
    
    timeIndex = str2num(datestr(dates,10))+str2num(datestr(dates,5))/12;
    plot(timeIndex,setup.macros(:,i),'-k')
    hold off
    axis tight
    %set(gca,'XTick',setup.timeIndex)
    if i == 1
        title(['Inflation']);   
    elseif i == 2
        title(['Output gap']);   
    end
end

for i=1:setup.nn+setup.nm
    subplot(numRow,numCol,idx+i)    
    hold on
    timeIndex = str2num(datestr(dates,10))+str2num(datestr(dates,5))/12;
    seriePlot = 12*output.factors(i,:);
    minY      = min(seriePlot)*0;
    maxY      = max(seriePlot)*1.2;
    X  = [timeIndex(1,1):1/12:timeIndex(end,1),fliplr(timeIndex(1,1):1/12:timeIndex(end,1))];
    % For shading above zero
    Y1 = [maxY*regime_2' 0*regime_2'];
    fill(X,Y1,[0.8 0.8 0.8]);
    % For shading below zero
    Y2 = [minY*regime_2' 0*regime_2'];
    fill(X,Y2,[0.8 0.8 0.8]);
    
    
    timeIndex = str2num(datestr(dates,10))+str2num(datestr(dates,5))/12;
    plot(timeIndex,12*output.factors(i,:),'-k')
    plot(timeIndex,12*(output.factors(i,:)'-1.96*sqrt(squeeze(AVarTheta1.VarFactorsHomo(i,i,:)))),'-r')
    plot(timeIndex,12*(output.factors(i,:)'+1.96*sqrt(squeeze(AVarTheta1.VarFactorsHomo(i,i,:)))),'-r')
    if i == 1
        plot(timeIndex,12*output.factors(i,:)*0+1,'--k')
    end
    hold off
    axis tight
    %set(gca,'XTick',setup.timeIndex)
    title(['Factor ',num2str(i)]);   
end

%% Plotting the factor loadings
figure('Name','Factor loadings','NumberTitle','off');
model    = output.model;
ny       = model.ny;
x2       = output.factors;
macros   = setup.macros';
factors  = [macros;x2];
[nx,T]   = size(factors);
jacobian = zeros(T,ny,nx);
for t=1:T
    x = factors(:,t);
    jacobian(t,:,:)  = model.gx + 2*model.gxx*kron(eye(nx),x);
end
jacobian_regime1 = prctile(jacobian(regime_1==1,:,:),[5 50 95],1);
jacobian_regime2 = prctile(jacobian(regime_2==1,:,:),[5 50 95],1);
subplot(2,nx,i) 
idx = 0;
names = {'Inflation','Output gap','weight inflation','weight output gap'};
for i=1:nx
    idx = idx + 1;
    minY = min(min(jacobian_regime1(1,:,i),jacobian_regime2(1,:,i)));
    maxY = max(max(jacobian_regime1(3,:,i),jacobian_regime2(3,:,i)));
    subplot(nx,2,idx) 
    hold on
    plot(setup.matSelect/12,squeeze(jacobian_regime1(1,:,i)),'--b')
    plot(setup.matSelect/12,squeeze(jacobian_regime1(2,:,i)),'-b')
    plot(setup.matSelect/12,squeeze(jacobian_regime1(3,:,i)),'--b')
    hold off
    axis(([setup.matSelect(1)/12 setup.matSelect(end)/12 minY maxY]));
    title(['Expansion: ',names{i}])
    
    idx = idx + 1;
    subplot(nx,2,idx) 
    hold on
    plot(setup.matSelect/12,squeeze(jacobian_regime2(1,:,i)),'--r')
    plot(setup.matSelect/12,squeeze(jacobian_regime2(2,:,i)),'-r')
    plot(setup.matSelect/12,squeeze(jacobian_regime2(3,:,i)),'--r')
    hold off
    axis(([setup.matSelect(1)/12 setup.matSelect(end)/12 minY maxY]))
    title(['Recession: ',names{i}])
end


end